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ABSTRACT 

Nuclear magnetic resonance is widely used as a probe 
of pore geometry and fluid composition in well logging. 
One of the critical assumptions often made is that the dif- 
fusion of fluid molecules is sufficiently fast to warrant 
the condition for direct mapping between the surface- 
enhanced relaxation rate and the pore geometry. In pores 
satisfying such a condition, but having a significant spa- 
tial variation of surface relaxivity (p), one can show that 
the one-to-one mapping may break down. The degree 
to which the NMR logging interpretation is affected has 
not been systematically studied until now. Extra relax- 
ation due to diffusion in internal field gradient is another 
example where spatially varying relaxation strength may 
obscure the direct relationship. Better understanding of 
their interplay may be exploited to our advantage. In this 
work, we theoretically investigate the interplay between 
the pore geometry, internal field and the inhomogeneous 
surface relaxivity. We develop a perturbative framework 
and compare its results with exact solutions obtained for 
a class of p textures in a pore with simple geometry. Its 
effect is quantified for a wide range of diffusivity and/or 
p strength. The result allows us to set the bounds for the 
change in the final slope of the relaxation curve and may 
serve as a useful guide for logging applications in real 
rocks with a wide range of pore sizes and fluid diffusiv- 
ity. We further employ large scale numerical simulations 
to perform virtual experiments on more complex situa- 
tions. Internal field and its gradient distributions were 
obtained and analyzed for up to 1.5^cm^ based on tomo- 
grams of carbonate rocks. We find that the texture of p 
based on the internal field gradient induces a small, but 
observable shift, compared to a random noise generated 
texture for which no shift is observed. 



PORE GEOMETRY, dp AND NMR LOGGING 

There exist potential pitfalls in the way NMR logs are 
interpreted ( |Kleinberg,1996| ). While it is widely agreed 
that the method is robust for simple types of porous me- 
dia, key assumptions for its successful application may 
become compromised progressively as their geometrical 
and lithological properties become complex. To be pre- 
cise, there are three necessary assumptions for the sim- 
ple mapping between the so-called T2- and the pore size- 
distributions to work: (1) The pores are practically iso- 
lated or periodic so that diffusive coupling( C ohen, 19821 
I de Gennes,1982l |Zielinski,20()2l ) among the pores may 
be neglected. (2) Within each pore, the so-called /a^r dif- 
fusion criterion is satisfied so that the relaxation is con- 
trolled by the weak surface relaxation strength, which 
will be represented as p (po, if uniform) from now on, 
rather than by the diffusive flux. The condition is given 
in terms of the control parameter sls = poL/D <C 
1 assuming a cubic pore of volume V = and dif- 
fusivity D of the fluid. ( [Brownstein, 1979| ) (3) The map- 
ping is based on the assumption, often made without any 
quantitative justification for a given rock, that p is uni- 
form across the pore-grain interface. Extensive investiga- 
tions were made on the first and the second( McCall, 1 99T\ 
Wilkinson, 1991! "Bergman, 1995! "Ryu,2001| |Ryu,2009a| 
Zielinski,2002,.Grebenkov,2007 ) issues, but the third has 
received relatively scant attention( Ryu,2d08{ |Ryu,20Q9b{ 
Ryu,2QQ9a| |Ams,2QQ61 |Valfous kaya,2QQ6 ). Unless one 



incorporates all these issues on an equal footing, it be- 
comes difficult to gauge uncertainty in an NMR log in- 
terpretation. 

The aim of this paper is to investigate systematically 
the nature of their violation and seek quantitative bounds 
for their experimental signatures, should they occur. This 
is done by first considering a simple pore with a class of 
p{r) textures which allows exact solutions for nontriv- 
ial cases. For a realistic pore geometry, we derive the 
pore geometry from 3D tomo grams(,Sheppard,20Q 4 ) and 
run random- walk simulations ( Ryu,2Q09b| undeF a set of 
controlled spatial profiles of p. Our focus is on gaining 
better understanding of whether and how the intertwined 
p texture and the pore geometry would make the issues 
acute or negligible. Careful experimental characteriza- 
tion of p is an invaluable component toward ultimately 
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increasing the utility of NMR logging for challenging en- 
vironments and such effort is emerging. ( [Keating, 2007 ) 
Our theory makes quantitative analyses possible in such 
investigations as well as guiding log interpretation when 
such information is unavailable. 

The evolution of the polarized proton spin density may 
be analyzed effectively using the eigenmode analysis of 
the underlying Helmholtz problem. ( Brow nstein, 1 979[ 
Ryu,200l| |Grebenk ov,2007|) For continuity, we will em- 
ploy the notational convention used in our recent paper. 



( Ryu,2009a ) Within the framework, the evolution of the 
total polarization M{t) is viewed as a superposition of 
the eigenmodes {^p(r)},(p = 0,l,...) each with relax- 
ation rate Xp. For <C 1, the slowest mode with p = 
dominates; viewing a porous medium as an ensemble of 
pores with a distribution P(Ao) of Aq values, one arrives 
at the backbone of current NMR log interpretation. 

Even though validity of the first two assumptions de- 
pends critically on the strength of po, its value, however, 
is often unknown. A popular method to estimate the 
po strength for a given porous medium is based on the 
assumption that the magnetization decays exponentially 
with its rate Aq ~ poS/V. From an NMR probe-derived 
T2 (~ 1/Ao) distribution, one may obtain an average 
< Ao >, and combined with the < 5'/^ > value (inde- 



pendently measured such as from BET( |Kleinberg, 1 999 



Keating,2007)), the estimated po strength, po, may be ob- 
tained: Po ^< Ao > / < S/V >. To illustrate a po- 
tential pitfall in this approach, consider Figure [T] which 
shows how actual Ao varies as po (and therefore hz) in- 
creases in a simple cubic pore of size V = L^. While the 
linear relationship Ao oc po holds for <C 1 as indicated 
by the borken line, the curve bends with an upper bound 
< Aoo = Djj. This means that the apparent strength of 
Po as estimated by the above method, will then be upper- 
bounded. 



P^<D 



6L* 



(1) 



For L ~ 100/im, and D = 2500/im^/sec, this means 
Po ^ 40/im/sec. As such an apparent po value is bounded, 
the apparent k, parameter 



'D 



6 



(2) 



is also upper-bounded, which gives the erroneous impres- 
sion that the premise < 1 is never grossly violated. Us- 
ing this circular reasoning, one would hardly encounter 
an empirical p value that would suggest the second con- 
dition is potentially violated. More accurate estimation 
of Po may be obtained if one uses the initial slope of the 
time-domain relaxation, which satisfies 



1 



po 



V 



(3) 
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Figure 1 : Behavior of the slowest relaxation rate Ao as a func- 
tion of the control parameter n = po^/^ in a cubic pore of 
V — with uniform po . The inset to the first panel shows the 
effective pore size parameter i^^jL. The value of Iq/L = 1/6 
corresponds to the inverse surface-to-volume ratio of the pore. 
The bottom panel shows the relative fractional contributions to 
Ao from the surface-integral (^Ao,(p)) and the volume-integral 
(^Ao,(i:))) contrbutions. (Eq|4]) 



for all values of k.. However, the range over which it 
holds become extremely narrow as n increases, and there- 
fore impractical for typical borehole applications. That 
leaves the rock-specific, accurate evaluation of po still an 
unresolved issue and we should be mindful of the poten- 
tial lapse in i\\tfast diffusion condition when k > 0.1 is 
observed. 

The extended nature of the pore space in natural media 
brings further issues. Suppose we manage to decompose 
such pore space into effectively isolated pore elements. 
How can one judge the robustness of the second assump- 
tion given these elements when their range of shapes and 
sizes vary significantly? For this, we first need to decide 
how to rephrase the second condition for general pore 
shape. The slowest relaxation rate Ao is composed of two 
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components involving the slowest mode (j)o (r) ,( |Ryu,2Q09a ) 

j dcj(i)l{v)p{v)+D j |V0o(r)|^(ir = aAo,(p)+aAo,(D) 

(4) 

and we note that their relative fraction controls the degree 
to which the second condition is met. When 9Ao,(p) ^ 
9Ao,(D), the relaxation is limited by the weak po, and 
is dominated by the slowest mode alone, validating the 
T2— pore size mapping. In the opposite limit, the slow- 
est mode profile acquires significant spatial variation, and 
higher modes gains weight, leading to a multi- exponential 
characteristic (B rownstein,1979 ) which translates into a 
significant spread in the T2 distribution. The bottom panel 
of Figure [T] shows how they behave over six decades in 
n{= pqL/D) for a cubic pore. The crossover takes place 
over the range near ~ 5. For pores with a shape that 
defies characterization using a single length scale, the 
length scale defined as 

may be used for the /^^ = p^l^/D parameter. It depends 
on the spatial profile of the slowest eigenmode and for 
cubic shape, it varies between L/6 (i.e. V/ S) for ~ 
and L/Afor n 00 as shown in the inset to the first 
panel of Figure [T] In general, its asymptotic values in 
both n ^ {) and 00 limits should be characteristic of 
the pore shape. Using this definition, we can put Aq into 

Ao = SAo,(p)(l + po4/l^), (6) 

an expression we find useful later (Eq. [58] ). 

For a collection of isolated pores with varying shape, 
one would proceed to apply the fast diffusion criterion 
via pqIq/D <C 1 for each individual pore. Isolated, 
large vugs, if present, may obviously incur exceptions 
to the direct mapping. There are subtle effects beyond 
this in realistic porous media. Significant spatial fluc- 
tuations present in the local diffusive coupling strength 
( |Zielinski,2QQ2| ) lead to enhanced inhomogeneity in 0o (r) . 
This means that 9Ao,(d) of Eq|6] gains more weight, ^0 
gets bigger, and the multi-exponential characteristics of 
the relaxation more pronounced. As a result, the sim- 
ple relationship between Aq and local pore geometry be- 
comes obscure. When the first and the second condi- 
tions break down simultaneously, the combined effect is 
of pure geometrical nature, and it justifiably invites active 
current research on issues such as NMR log-prediction of 
permeability. 

Unfortunately, the difficulty may not stop at the purely 
geometrical level. Suppose we have a system with a uni- 
form Po satisfying the first two criteria, for a range of 
pore sizes a G [amin-,cimax\' Then the observed rate 
distribution Prate (^o) leads directly to the size distribu- 
tion Psize{o)- It is hypothetically possible, however, to 



have a collection of isolated pores of a suitably chosen 
size ao, with a distribution Pp{po) of po values assigned 
to each, that will yield the identical relaxation behavior. 
Since there exists scanty empirical data on p(r), and it is 
not always possible to have the actual 3D pore geometry 
known, it is not clear how to quantify uncertainty orig- 
inating from break down of these assumptions respec- 
tively. Accuracy of any physical property dependent on 
Ao and its probability distribution, such as permeability 
and porosity, will then be controlled by the uncertainty 
present in Ao. 

To bring some clarification to this, we recently devel- 
oped a theoretical method to solve for the change in the 
relaxation rates {Xp} and their associated eigenmodes 
{^p} which sets the bound for uncertainty for a prop- 
erty derived from Ap's.( |Ryu,2008] |Ryu,2Q09a| ) In sum- 
mary, we consider the spatially fluctuating part of the 
p{r) given by 

6p{r) = p{r) - Po (7) 

where po is the average of p{r) over the interface, and 
derive the fractional shift in the slowest rate taking Sp/po 
as the perturbation, 

in terms of various surface overlap integrals that involve 
the eigenmodes of the uniform po case, 

5p,o^ = j da4>l{v)5p{v)4>l{v). (9) 

In the case of a discrete hemispherical (5p,( |Ryu,2QQ9a| ) an 
exact solution was also obtained, extending the bound for 
cases where 8p/ pQ is not small. Due to high symmetry 
of the spherical pore, it left room for other simple cases 
where deformations in pore geometry and 5p patterns are 
more explicitly intertwined. 

In the following, we therefore extend the results to a 
rectangular pore. A set of non-trivial ^p's, appropriate for 
setting the bounds, are considered. For a rule-of-thumb 
type bound for log interpretation, we lay out steps to es- 
timate the degree of change in 5\q/\q (i.e. 5T2/T2) ex- 
pected for a combined geometrical and p— textural devia- 
tions from a pristine condition assuming that the porosity 
is preserved. We also demonstrate by numerical simu- 
lations the importance of relationship between the 5p{y) 
texture and the ^o(r) profile and their symmetry. The 
manifestation of such effect is first considered for quadra- 
ture patterns of (5p in a cubic pore followed by cases with 
tomogram-derived 3D pore geometry. In the latter, di dp 
texture based on the internal field and another based on 
the correlated random noise sequence were imposed on 
the interface of a carbonate pore matrix. 
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RECTANGULAR PORE : BOUNDS FOR p(R) AND 
GEOMETRICAL VARIATIONS 

Defining pore geometry and p textures 

Once we allow general variations in (5p as well as pore 
geometry, the phase space quickly expands beyond our 
means. To make the problem tractable, yet meaningful, 
we consider a situation in which both vary under con- 
straint in a space of dimension de . Consider a rectangular 
pore with the porosity function: 



1 ifrc, G [-^, ^] for alia = l,...^ 
otherwise 



(10) 

with its geometry tuned by aspect ratios La/Lp preserv- 
ing the total volume V = YlLa- Such a variation fails 
to incorporate the heterogeneity present in natural me- 
dia, but it is a good starting point. We systematiclly ex- 
plore how 5p interferes with the pore geometry, specifi- 
cally a porosity-preserving distortion and eventually a di- 
mensional crossover, and to what degree they affect the 
T2— distribution. 

For the p texture, we consider a situation where each 
side of the rectangle may have a distinct value pcx± where 
a± indicates the left (— ) and the right (+) interface in the 
a— direction . For each dimension, we introduce their 
average for each direction: 



Pa 



Pa+ + Po 



and the asymmetric variance 

_ Pa-\- — Pa- 



Spo 



(11) 



(12) 



Generally, pa varies for each direction a. Define gross 
average of p^s across the whole interface of the pore as 



<p>-- 



(13) 



where Sa = n/3^a Lpis the cross-sectional area normal 
to fia . It is more convenient to use dimensionless control 
parameters, ctq, and e^'s. defined as 

Pa-<P> .... 

—rr- (14) 

< p> 

is the relative deviation of each planar value Pq, from the 
average, may vary for each direction, but they consti- 
tute symmetric part of the p fluctuation as Pq,± 's are equal 
in opposing planes unless cTq, 7^ 0, where 

5pa 

Pa 

represents the asymmetric fluctuation of the p between 
opposing planes. Note that 



(15) 



(16) 



by definition. Note that and cr^ thus defined are analo- 
gous to the compression/dialational and the shear defor- 
mations in the theory of elastic deformation. Finally, we 
define 

_ PaLa 

l^a = (17) 

as an analog in each dimension to the n parameter. 
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Figure 2: Schematics of three cares of a rectangular pore. For 
clarity, 2-dimensional cross-sections are shown. The strength 
of the gray scale on the pore wall indicates the strength of p(r) 
on each plane. Left panel shows the symmetric rectangle (center 
panel) with (ja = (for all a = x, 2/, ^); the center shows an 
asymmetric rectangle with cr^ / 0. The right panel shows the 
base system, in which all p^'s are equal io pg —< p > and 
Ca, = in all directions. The volume of the base cube is 
chosen to match that of the rectangle. 

We obtain analytic solutions for the slowest relaxation 
mode and its eigenvalue Aq for each of the three cat- 
egories schematically described in Figure [2] The case 
with ctq; = for all a will be referred to as the sym- 
metric texture and will be indicated by super/subscript s\ 
Cases with non-zero cFoc values, (but with may or may 
not equal 0) will be called the asymmetric textured with 
super/subscript a where appropriate. The globally uni- 
form case is a special instance of the symmetric texture 
for which poc are all equal and will be indicated by su- 
per/subscript g . With the symmetric texture, there is no 
finite minimum separation between a pair of interfacial 
points with distinct values of p's. Due to the separability 
of the coordinates, the solution for such situation can be 
trivially constructed out of the one-dimensional problem 
with the same boundary condition at both ends. We also 
impose 

Pc. > 0; 1 > cTc, > 0; (18) 

without losing generality. For symmetric textures, it is 
further required that (7^ = for all a. The globally uni- 
form case requires also that Cq, = for all a. ctq, = 1 is 
rather special in that it forces the p^- to vanish on one 
of the opposing planes with an interesting consequence, 
(see Figure |7] and the discussion following it) 

To incorporate variations finer than described above, 
one may divide the rectangle into a set of smaller ones. 
The minumum size for these sub-rectangles will be set 
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by the smallest length scale dictated by that of Sp{r). 
Each of these sub-rectangles couples with its adjacent 
neighbors through the diffusive coupling, leading to a 
large set of coupled equations instead of a simple 2 by 
2 homogeneous equation that we have here. One may 
also consider a sphere (or a cylinder) or ellipsoidal ge- 
ometry. We have worked out a perturbative scheme for 
the spherical geometry (Ryu, 2009 a) and obtained exact 
solution for a simple binary hemispherical texture of p 
( |Ryu,2009b) ). Extension for more complex angular tex- 
ture of is straightforward in theory, but involves a large 
number of angular modes. Rectangular (or cube) geom- 
etry makes the numerics much simpler for the class of 
p textures of Figure |2] while at the same time allowing 
a nontrivial geometrical variation. Note that the sym- 
metric texture for a rectangle (or cube) corresponds to 
the L = 2m (m = 1,2,...) spherical harmonics modes, 
while the asymmetric texture would bring in odd harmon- 
ics, L = 2m + 1, of the spherical pore. 

General Solution 

The boundary condition for the general p texture defined 
in Figure |2] now leads to the pair of conditions in each 
dimension a: 

(L,/2) + D\/aMLa/2) = (19) 

Pa-^0,a(-^a/2) - DV a(t>0 ,a{- L a / 2) = 

for the lowest eigenmode (/)o({rQ,}) = </>o,a(^a) with 
each (/)o,q; in the following form: 

^o,a(^a) = cos(A:c,rc,) + 5^ sin(A:c,rc,) (20) 
with the normalization condition 



/ drocCos^ {kocra) + I dr a cos^ (kaV a) = I 

(21) 

The constrained Sp textures and the geometry render the 
problem separable in each dimension, so we have a one- 
dimensional problem for each a satisfying: 



JC 







(22) 



with the matrix JCa given by 



cos(^) + sin(^) 



V cos ^ - ^qo. sin ^ a, sin ^ 

and the dimensionless wavevector qa 

Qa = kaL^/ TT. 



(23) 



(24) 



For a solution to exist, the determinant of JCa should van- 
ish, therefore the solutions are given by the null space of 



ICa . The solution space is spanned by eigenmodes cor- 
responding to an infinite set qa{i)j ^ = 0, 1, 2, . . . each 
satisfying 

— (2cos — l)+(l-( ) -cr^)cos— -sm— - = 0. 

(25) 

The smallest of such {qa}, Qa{^) in each dimension con- 
tributes to the rate of the slowest mode as a function of 
sets of independent parameters z^:, a, L in each dimension, 

Aa,0(W}, W}, {L4) = V Aa,a(0) = Vi^(^)'^^(0) 

Oi ex. 

(26) 

and takes the form of a set of parallel channels. Introduc- 
ing the diffusion time 



it can be put into the following form: 



Solutions with a Symmetric Texture 
For a symmetric texture ({ctq, = 0}), /Cq, is 
/ 



(27) 



(28) 



-qa sin ^ 



qo, cos ^ + sm ^ 




(29) 

and the boundary condition factorizes. For the slowest 
rate, the symmetric solution should be taken 



h{{ra}) = W^oi cos{kara) 



with the normalization condition 

/ draCos^{kara) = 1. 
The boundary condition yields 

TT 



(30) 



(31) 



cos((7c,7r/2) 



-qaSm{qa7r/2) = 0. (32) 



This leads to the class of solutions equivalent to those 
used by Brownstein and Tarr for the simple geometry. 
( iBrownstein, 1 979 J The rate for the slowest mode is then 

A.,o({«4,K}) = E— ■ (33) 

a 

Note that it takes the form of three competing diffusion 
channels each of rate ^ weighted by q'^ ^ factor. It is 
useful to further examine limiting behaviors in this sym- 
metric case: In the limit of K,a ^0, 



\im qa{0) =q^ 



/2k.. 



TT 4 

(34) 
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while in the opposite Hmit /^^ ^ oo, 



lim (7c, (0) = q^ = 1 



(35) 



where > 10 is required for the expansion in ^ to be 
accurate. In the special limit where all a^q, <c; 1, 



lim Xs,o{{kcx}, {ra}) ^ ^^(1 



(36) 



Note that this approaches (< p > — < pn > 
Under this condition, this form shows that the symmetric 
texture leads to that of the uniform case. In the opposite 
case (WhZa ^ 1), 

1 - 



lim Xs,o({kc}, {r^}) 



For a general symmetric case where either of the limits 
cannot be taken for all directions, i.e. when the hz^ Pa- 
rameters fall in the intermediate zone, appropriate inter- 
polations may be made. 




1 
logio/^ 



Figure 3: Qp as a function of k, and a. The inset show a 
3D plot of Qp{K,,a). Note the rapid drop from Qp ^ 1 to 
0.5 toward the large k and a ^ 1 corner, as it corresponds to 
doubling of the wavelength for the slowest mode and accom- 
panying phase shift. Upper panel shows the series of curves 
logio Qp loSio ^ values of a = 0., 0.1, 0.2, ... 1.0. In 
the main panel, all curves converge to sl Qp ^ k behavior for 

K < 0.1. 



Cube with Uniform p 

The solutions obtained above for general Ca-cTa — La 
variations, may be compared against that of a pore with 
globally uniform p = pg and cubic geometry with the 
size, Lg"" that is equal to that of the rectangle Yl^"^ La- 
Defining the effective Kg 

Kg=<P> (38) 



it is related io via 



2V 



D 



f^a VLg 



^Ll S 



(39) 



using the fact that SaLa = V for any a for a rectangle. 
The slowest rate for the globally uniform cubic pore takes 
the simple form: 



(40) 



where g^(0) is the smallest of the g^'s that satisfy: 



cos((7^7r/2) Qg sin(^^7r/2) = (41) 



and we introduced 



(42) 



(37) Instead of the volume matching criterion, had we chosen 
Lg = 2de{S/V)~^ so that the cube has the same surface- 



to-volume ratio to that of the rectangle, we recover the 
expected behavior that Xg^o{K,g) Aq (a^q, , 0) of the sym- 
metric case when hiajf^g <C 1 and pa =< p >• How- 
ever, in many contexts, it makes more sense to impose 
the equivalent- volume criterion, i.e. deformation of pore 
geometry while preserving overall porosity. In the fol- 
lowing, we will employ such a criterion with = V. 



Recipe for Estimating Bounds 

In summary, we note the following: For a cube of side 
length Lg with a uniform pg, the rate, A^^o, is controlled 
by the one-dimensional condition (Eq|4T]) weighted by 
the arithmetic mean of the geometrical factors i.e. diffu- 
sion rates (Eq|4Q|). In the symmetric case (rectangle with 
Ca's), the boundary conditions become distinct for each 
direction (Eq|32|), controlled by the hZa factor; the slowest 
rate, Xs^o, is given by Eq[33|and is a function of K,a and 



Tq; 's. The uniform pg is related to the pa 's through Eq[T3| 
In the asymmetric case (ctq, 7^ for at least one dimen- 
sion), the even and odd modes mix (i.e. phase shifts) with 



its degree controlled by ctq,, Eq 25] The resulting rate for 
the slowest decay mode (Eq|28), Aa,o, is a function of 



well as cTq^'s. The diffusive rates r^, ^'s control 
the purely geometrical aspects of the pore. La • 

If one is interested in the change of the rate in going 
from the base (uniform, isotropic) cube to the rectan- 
gular pore with a symmetric texture, its fractional shift. 



5X 



g:s,0 



X 



A. = 



^,0 
5X, 



g:s,0 



As,o,is given by 



(0). 



(43) 



where the subscripted Qs^a is used to indicate that it is the 
solution of the symmetric case. For comparison between 
the asymmetric and the reference cases with SXg:a,o = 
Xg^o — Aa,o, we obtain 



A, 



SX 



^9,0 



1 



Ed 



)■ 



(44) 
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Figure 4: as a function of k. It is vanishingly small for 
K, < 0.01 but increases toward the value of 1 for large k. 

These expressions treat the impact of both geometrical 
and p-textural changes on an equal footing. Separation 
of their individual impact is straightforward: For geo- 
metrical impact alone, one may take the first expression 
evaluated with = for all a for a given rectangular 
geometry (let us identify this special case with subscript 
s with Pa = pg for all a, and = PgLa/D). Compar- 
ing this to the cubic, uniform p case. 



A. 



,0 4 Ta qm ^ ^ 



For the difference between the p-textured and a uniform 
p with the same rectangular geometry, use the 5— case 
as the reference instead; the fractional shift between the 
uniform and textured rectangle is 



A, 



(1 



(46) 



Due to the weighting factor )^s,a{^) / '^s,a{^)^ the 
largest As,q,(0) will dominate, which in the s— system, 
is equivalent to the shortest La since pa =< p > for 
all a. Therefore, for an extremely anisotropic geome- 
try, Li <C La {a = 2, . . . de), such as a slab-like pore, 
this reduces to (1 — q'^ i{^)/Q'i dominated by the 

p— variation in the most constricted dimension. Note that 
while ^A^:a,o = ^^g:s,o + ^^s:a,o IS truc, onc cannot er- 
roneously assume A^ = A^eom + Ap. 

For a general variation of {pa, cr^, Lq,} (a = 1, . . . dg) 
values, we have a 3 x dimensional phase space. How- 
ever, one can determine the impact of moving in such a 
space by reading a few numbers off a universal function 
Qp{a^K,). First, pure geometrical aspects are incorpo- 
rated in terms of rQ,'s once the aspect ratio of the pore is 
set. Next, more subtle aspect involving both the geome- 
try and the p texture is addressed via the series of master 
curves Qp{i^, cr) from the smallest \Q\ value that satisfies 



the condition 

-(2cos^^-l) + (l-(^f- 



7T 



2. Qtt . Qtt 
- a ) cos — — sm — — = 0. 
^ 2 2 

(47) 

This establishes a manifold in the z^:— a — Q space, plotted 
in Figure [3] that determines the slowest rate for any con- 
figuration prescribed in Figure |2] Obviously, one can fur- 
ther construct manifolds that correspond to faster modes, 
although logistics of tracking among closely spaced higher 
eigenvalues may not be trivial in practice. Both qg{0) 
and qs,a{^) of Eq 43 can be read off from such Qp at re- 
spective values of K, = K,g(Eq\S^ and = /^^^(EqjTT]) with 
(7 = 0. Likewise, the values of qa,a{^) are given from 
Qp at K = hZa(^^^]} and a = ctq, (Eq[T5]). 




Figure 5: Top: Ae as a function of n. It is non-negligible 
for small k values (^ —1) and becomes negligible for large 
K > 100. Even for small however, the overall contribution 
from Ae averages out due to the condition ~ ^' ^^^^ 

that the functional form is also essentially identical to A^, re- 
flecting the fact that the symmetric p texture (i.e. Ca) has an 
effect similar to that from geometrical dilation/contraction of 
the pore dimension. Bottom: A^^ as a function of k and for val- 
ues of cr = 0.01., 0.02, 0.04, 0.06, 0.08, 0.1, 0.2, 0.3, 0.4, 0.5. 
The inset shows A^/a. For a < 0.1, the curves tend to con- 
verge toward the universal form for the whole range of k. 

As an application of this recipe, let us consider three 
hypothetical variations: one in geometry (tq,) alone and 
the others in terms of Ca and cTq^'s in comparable frac- 
tions, say 50% to get an idea of their relative signifi- 
cance. (I) For a pure pancake-like geometrical deforma- 
tion specified by 2r^ and Ty^Tz rg/\^, the 
rate becomes faster with A^eom = —0.222 {hZg = 0.32), 
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-0.286 {i^g = 1.6), -0.435 {f^g = 6.4). Now, con- 
sider a cubic pore with variations in p only: (II) For a 
symmetric Sp texture with Cx = —0.5, Cy^e^ = 0.25, 
the rate slows down with = 0.0062(/^^ = 0.32), 
0.024(a^^ = 1.6) and 0.039(a^^ = 6.4). These small val- 
ues are due to cancellation of larger numbers from each 
direction (e.g. with k, = 6.4, we had A^^^; = 0.1016 and 
As^y^As^z = —0.0311). (Ill) For an asymmetric texture 
with a a = 0.5 for all a's, it always slows down with 
Aa = 0m4{i^g = 0.32), 0m9{i^g = 1.6), 0.106{i^g = 
6.4). 

Small K limit 

Without going into specific excursions in the 3 x dg- 
dimensional phase space, we can make general observa- 
tions on how changes in pore aspect ratio and p textures 
affect the rates. Consider a rectangle of {Lex} {S/V = 
J]]^ 2/I/c^) with an asymmetric textured {pa} and {(Jo} 
values. Consider also its variant, the symmetric system 
in which VcTq, = 0. We also construct an equivalent 
cube of side Lg and globally uniform < p >. From 
Eq|40| the rate for the globally uniform cube is deter- 
mined by Qp{hZg^O). Introduce aspect ratio variations, 
then the fractional difference in the rate (Eq|43]) between 
the symmetric and the cube is controlled by the (1 — 
y-Ql{^a, 0)/Q l{i^Q , 0)) factors. In the limit where all 

\2 l+'^a/iA 



A>:'s are small (Eq 



and using Eq 38 



34), they become (l-(f5)" 
we can show 



l+/^«/4 



lim A 5 



1 

dp. 



-(l + e«) + 0(A>:,). (48) 



Furthermore, due to the property Eq[T6j the term linear 
in to. vanishes. Therefore, we have 



lim As 



1 



d. 



-y 



9 

La 



0{Kg). (49) 



where the last equality follows from the porosity pre- 
serving definition Lg = V^l^-. The expression above 
(and more generally Eq|43]) reveals the separation of the 
geometry (through factors) and the class of ^p(r) 
variation (via vanishing of the e^/Lo, factors), at least 
in the small Kg limit under the fixed volume condition. 
Note that with e^-type variations, the slowest modes are 
still all symmetric functions in each dimension. Volume 
preservation requires that both compression {ta > 0) and 
dialation(eQ; < 0) should be present in different direc- 
tions, and as a result, their impact cancels out to first or- 
der. Even if it survives, it is eclipsed by the impact of 
change in the surface-to- volume ratio. For large values 
of Kg, terms involving e^s may survive to contribute, 
and As should be numerically evaluated out of the mas- 
ter curve Qp{K^O) evaluated at k = Kg and a^^'s, each 
multiplied by the geometrical factor, (^)^. 



For the asymmetric texture (a a / 0), the mode pro- 
file, in addition to the contraction/dialation, goes through 
a phase-shift in general. (Ba 7^ in Eq|2Q|) as it is being 
pulled at one end and pushed at the other by the asymmet- 
ric pa±- (See also Figure [7] and discussion below) The 
fractional change in the rate is constructed from read- 
ing off Qp curves, and it involves moving across varying 
(Jc^'s, i.e. excursions in the vertical direction in Figure[3] 



!2^[^(l+e),0] 
1 0.0 




Figure 6: Top: Fractional shift factor for the symmetric texture 
with (ja = and — 1 < e < 1. By definition of < p > and 
Ca, we have e = for 1 -dimensional system, and for de > 2, 
^^(ca/La — SLa/L"^) = should be satisfied for a general 
symmetric case. Bottom: Fractional shift factor for the texture 
with 1 > (Ja > 0, = 0. No assumptions are made about the 
size of a. 



Large k limit 

Let us now examine geometrical and textural deforma- 
tions off an arbitrary value of Kg which is not necessarily 
assumed small. For small fractional changes in k^^ aris- 
ing from a geometrical deformation alone, Eq|45]yields 



A 



geom 



1 X - dL, 



-(1 + AlM). 



where 



Kl{k) = 1-2k 



dQp{K^ 0)/dK 



(50) 



(51) 



Qp{k,0) 

which is plotted in Figure |4] In the small k limit, it 
converges to the value of 0, making Ageom converge to 
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Figure 7: Changes in the profile of the slowest eigenmode 
00 (r) for different values of k and a. We assume all and 
aa's are zero except in the x-direction for which ax / as 
indicated for the solid curves. The broken curves are for the 
uniform case with ax = 0. The shades indicate whether the 
profile increases (darker) or decreases (lighter) when ax / 0. 

SS/S as expected, while in the opposite limit, the con- 
vergence is toward the value of 1, so that it gradually ac- 
commodates the slow diffusion asymptote Aq,(0) cx 

In the case of changes arising from p texture alone, 
we take the s— state as our reference with Qs^a^ solely 
determined from pg and {La}, and with <C 1, cr^ <C 
1, 



A, 



where we define 



Ae(/^) = -2f^d\nQp{f^,0)/dK 



(52) 



(53) 



and 



A^(/^,cr) = -2d\nQp{f^,a)/da. (54) 

Figure[5]show the functions Ae{hi) and Acr(A>:, a) respec- 
tively that can be used for small values of and so 
that Taylor expansion of Qa^a around the q^^a is valid. 
Note that when plotted as Acr/cr (Inset to the second panel 
of Figure [5]) the curves for small a values (< 0.1) con- 
verge toward a universal curve which peaks around ~ 




*%(0,0,L/2)^ 




Figure 8: Top-left panel: Profile of 0o(r) for n = 1.6 for a 
cubic pore. The colormap was chosen to emphasize the deple- 
tion of 00 near the eight corners of the cube. The rest shows the 
profile of ((/)o)^ in one of the interfacial planes normalized its 
value at the center of the plane, k — 0.32, 1.6 and 6.4. 

uniform p face centered 5p cornered 5p 



Figure 9: Three p textures considered for comparison. The 
first panel shows a cube with a uniform po • Two non-uniform p 
textures are shown. The mid-panel shows face-centered 3p 
which has an enhanced p in the center of each face. The right 
panel shows the cornered 5 p which has enhanced p in the eight 
corners of the cube. The total area of the red is chosen to be 
equal in both cases. 

3.05. Since the contribution of ctq, texture to Ap is given 
in cr Act , this implies that the contribution to the shift be- 
comes second order in cr. This is consistent with the com- 
parison made earlier between the exact and the perturba- 
tive solutions in the case of a spherical pore (Ryu,20Q9b| ) 
and with the symmetry requirement, as the shift should 
remain independent of the sign of cr. 

Figure|6]shows the factors that control 5\q for the sym- 
metric (e 7^ 0, cr = 0: top panel) and the asymmetric 
(e = 0, cr > 0: bottom panel) cases while the geometry 
is held fixed. In the symmetric case, the contributions for 
the positive and the negative e's display strong asymme- 
try for large n. Note that we had observed earlier that 
these tend to cancel out for small n'^ due to the condition 
Eq[T6l For a large it is no longer the case. The second 
panel largely duplicates what we had found for the spher- 
ical pore.( ,Ryu,2QQ9a^ The peculiar evolution of this fac- 



9 



SPWLA 50*^ Annual Logging Symposium, June 21-24, 2009 



0.1 







III- 
uniform p 






face centered p - 






cornered p 






exp(-t /x g ) 




K = 0.32 


1 1 1 



0.00 0.05 0.10 0.15 0.20 0.25 
t [sec] 



0.1 



- 3.0 


















■ 








0.001 0.01 


: K=1.6 





0.00 0.02 0.04 0.06 0.08 0.10 
t [sec] 



0.1 



0.001 



1 1 3.0 

- 

\ 2.0 

- 

xxXo.o 


:,„iv, ,11, 




: K = 6.4 


'^O^OOl ^ " '^O.Ol ^ 



10 20 30 40x10" 

t [sec] 



Figure 10: Comparison of NMR responses simulated for the 
three p textures of Figure [9] The first three panels (each with 
hi = 0.32, 1.6, 6.4) show M vs. t with uniform p (green), face- 
centered 6p (red) and cornered 6p (blue). The broken lines for 
short t indicate the predicted initial slope of the decay curve, 
~ e~^^'^^ for each n. The insets for n — 1.6,6.4 show the 
Laplace-inversion (T2 -distribution) results with the same color 
scheme. Also shown in the insets are the positions of the slow- 
est T2 from exact calculation indicated by the short vertical bar 
in yellow. 

tor from a peaky structure to a step-like shape as a ^ 1 
can be understood if we examine the way profile of the 
eigenmode evolves as a increases. Figure [8] contrasts the 
profile (\)%{x) of the slowest mode with a = 0, (shown 
with a broken curve in all panels) and 0o(^) for finite a 
values (solid curves). The changes between ^^^^ ^0 
are shaded gray. The left column is for = 2, right col- 
umn with K — For small k and a, (top-left), the effect 
is generally a moderate shift in phase. As a increases, the 
wavelength tends to shrink. This is most pronounced in 
the large k — a values (right-bottom panel). Note that in 
this limit, one has a large p on one side, and vanishing p 
on the other. Therefore, the system evolves from where 
the span L of the pore matches the half-wavelength of 
(broken curve) to a highly asymmetric profile (solid 
curve in the bottom-right panel) in which L equals the 



quarter- wavelength of (^o- Doubling of the length scale 
in (/)o leads to a decrease in Aq by a factor of 4, as indi- 
cated by convergence toward 0.75 as cr 00 in 
the bottom panel of Figure [6] 

QUADRATURE p TEXTURE ON A RECTANGU- 
LAR PORE 

The classes of texture considered in the previous section 
capture basic aspects of p and geometry in their entangle- 
ment with each other. However, it misses a subtle ingre- 
dient that may play an important role in media with non- 
trivial geometry. Recall how we showed that the first or- 
der Ca contribution to (5A5^o(0) vanishes. That arose from 
the constraint that pa is uniform in each plane. For more 
general situations, it is not necessarily so as the symme- 
try of Sp and the eigenmodes on the interface play an 
interesting role. For the leading oder contribution to the 
shift, we showed: 

6X0 = £ ^a(/)g(r)(5p(r)(/)g(r) + 0{6p^) (55) 

with Sp constrained to satisfy 

£ daSp{r) = 0. (56) 

In the case of a spherical pore, we had uniform across 
the interface, and therefore this first order contribution 
was shown to vanish. It was further shown that the higher 
order contribution always acts to slow down the decay, 
i.e. Aq > Aq. In a rectangular pore, it is no longer the 
case for a general Sp and we have a chance to observe 
the first order effect if Sp{r) variation is chosen to have a 
non-trivial overlap with respect to 0o(^)- is further ex- 
pected that SXq may become either positive or negative. 

Let us demonstrate that such a case is readily observ- 
able with a rectangular pore. With a uniform po, its slow- 
est mode is readily obtained: 

«=1 V 2 U + ko,o.L^ ) 

with each G [— Lq;/2, Lq,/2]. The first panel of Figjs] 
shows the relative strength of in ^ cubic pore for (= 
poLa/D) = 1.6. The colormap is chosen to emphasize 
the depletion of (/)q near the eight corners of the cube. 
Note that within each interfacial plane, 5*^ with its sur- 
face normal fia, ^0 ^ ^^^^^ maximum at the center 
of the plane. The next three panels show the profile of 
(0q)^, as it appears in Eq|55j in one of the interfacial 
planes for three values of a^: = 0.32, 1.6, and 6.4 normal- 
ized with respect to its value at center of the plane. Now 
we consider three Sp textures as described in Figure[9]im- 
posed a cube of size L^. Two cases, one with enhanced p 
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Figure 11: Top panel: Fractional surface-integral factor of 
(00 )^ over the red area on which p{r) is enhanced. The rel- 
ative weights of Ir and h controls whether the rate increases 
or decreases with respect to the reference case. Their degree is 
further controlled by the /p = (1 + poio /D)~^ factor which is 
a purely geometrical parameter and is shown in the inset. The 
second panel shows the curves of Figure[TO|with time scaled by 
the factors (1 + 5\q/\q) evaluated using values of fp and Ir 
above . Total of nine curves are shown, three each (following 
the same color convention for each 5p texture for the three hi 
values. The three curves collapse to a single curve for each k. 

strength in the square centered on each plane (designated 
as face-centered), another with enhancement in the eight 
corners of the interface (cornered) are compared with re- 
spect to the uniform po- For the uniform case, the rate is 
given by 



^0 



the fractional shift, Eq|55] becomes 

<5Ao _ §^da<t>°oir){pir) - po)<t>°oir 



^ Po 



I, 



1) 



h Po 



h 



1 



where /p = 1/(1 + poio/D) and 
L 

6 



Ir 



I (58) 

(59) 

fp 

(60) 



is the surface integral restricted to the red part (S^) of 
the interface in the Figure [9] and similarly for I^. The 
particular choice we made as depicted in the figure leads 
to the average po = Pbf + where p^^b) is the local p 
value on the red (blue)-part, which takes up 1 /4 of each 
plane. 



Figure 10 shows the simulated relaxation curves ob- 
tained using the random walk method(Ryu,2008 ) for the 
three cases in each panel. Three panels with progresively 
larger values are shown. In these calculations, we fur- 
ther chose to have pr = 2pb. (therefore po = \ph) The 
overall trend is that the slowest rate for the face centered 
pattern got faster, while the rate for the cornered texture 
got slower compared to the reference case with the uni- 
form Po- The trend becomes more pronouced for larger 
The same can be observed in the T2 -distributions of the 
same data shown in the insets, although the amount of 
change may not look unambiguous for those experienced 
with the inherent uncertainty in such a representation. 
The trend in the corresponding time-domain curves is 
free from the inversion-related issue and is real. The bro- 
ken curves attached to each time domain graph at early 
times indicate the predicted exp(— t/rg) behavior where 
Ts = poS/V. I/ts is also indicated in the inset as a short 
yellow bar. 

For these calculations, the first order contribution to 



the fractional shift of Eq 59 is shown to be 

SXo 



+ /6 



5 ' 



(61) 



Note that in the Figure 11 the values of Ir / {Ir ^ h) ^ 
0.25 for the centered texture, while it is < 0.25 for the 
cornered texture. Thus the first order contribution de- 
creases the rate for the cornered and does the opposite 
for the centered. Reading the values of 7^7^ and fp from 
Figure 1 1 we obtain, for the cornered texture, ^Aq / Aq = 
-0.014295, -0.0465, -0.0542 for n = 0.32, 1.6, 6.4 re- 
spectively. For the face-centered texture, the rate is pre- 
dicted to increase with ^ = 0.0149, 0.0548, 0.08069 
for hi = 0.32,1.6,6.4. This is in excellent agreement 
with the simulated results, as indicated in the last panel 



of Figure 1 1 where we show that the curves for different 
textures from Figure [T0| all collapse when plotted with 
respect to the time scaled by the (1 + SXq/Xq), using 
the fractional shift values found above. The agreement 
hardly leaves any room for the higher order perturbative 
contribution to make any significant addition. 

INTERNAL FIELD-LIKE p TEXTURE VS. PATCHED 
p ON A CARBONATE ROCK 

So far, we have not addressed whether there should exist 
any correlation between the spatial profile of Sp{t) and 
the underlying pore geometry in natural media. When it 
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Figure 12: Top panel: Internal field calculated to first order in 
Ax- The sample is a carbonate rock (packstone) of distorted 
rectangular shape with dimensions 1 x 1 x 1.3 cm^. The ap- 
plied field is along the long axis of the sample. The calculation 
was done for the entire tomogram volume 1.5 x 1.5 x 1.3 cm^ 
including the free space surrounding the rock, which is mirror- 
reflected and then periodically repeated. The colormap was 
adjusted for optimal contrast of the features. The part for the 
surrounding water is displayed with enhanced transparency to 
lessen obstruction of view. Part of the rock was cut out to re- 
veal its inside. Bottom panels: cross-sectional cutouts of the 
internal field. The left panel is from the same sample as above. 
The right panel is from a Berea sandstone with cross-sectional 
dimension of 5.2^mm^. 

does, it would depend on the geological history of the for- 
mation; Extrinsic factors such as the presence of strong 
magnetic minerals and their distribution may vary in a 
haphazard manner from one area to another, and one can- 
not expect to have a generic profile that covers them all. 
Systematic experimental analyses are emerging only re- 
cently with relevant details ( Keating, 2007). In previous 
sections, we considered a few artificially imposed tex- 
tures: In the first section, it was uniform in each inter- 
facial plane, but allowed to vary from plane to plane. 
In the preceding section, patterns were imposed delib- 
erately to control degrees of conmiensuration with the 



00 (i*). In our earlier work, (Ryu,2008 ) we employed 
grain-specific assignments in a random beads pack as 
well as correlated-random noise (patchy) pattern with vary- 
ing correlation lengths. In a spherical pore, both exact so- 




Figure 13: 512^ volume (1.46^mm^) of a high resolution to- 
mogram of a carbonate rock used in NMR simulations. Blue is 
the pore space, white represents the grains. The volume con- 
sists of mirror-reflected images of the seed volume of dimen- 
sion 256^ as indicated by the shaded cutout. The lower panel 
shows the pore-grain interface only. 



for hemispherical assignment. Ams et al ( |Ams,2006l ) 
used grain-assignment in numerical simulations on pores 
generated from tomograms. Using random assignment at 



lution and numerical simulations were obtained ( ,Ryu,2009a| ) 



a voxel-level, Valfouskaya et al ( Valfouskaya,2006| ) ob- 
served negligible effect. Most of these observations can 
be understood within our theoretical framework yet there 
are sources of 5p that defy an easy categorization. As an 
example, we consider a case where the surface-relaxation 
may vary in strong registry with a aspect of pore geom- 
etry, yet it is not clear a priori whether it would work in 
its favor or against. 

The most widely accepted mechanism for the micro- 
scopic origin of p is based on the engagement of a mi- 
grant proton spin with a surface-embedded paramagnetic 
ion spin. ( Korr inga,1962[ [Brown, 196H |Kleinberg,1996| ) 
On a larger length scale, the pore matrix and its filling 
fluid have different magnetic susceptibilities (let us de- 
note the difference as Ax), and this gives rise to an inter- 
nal field B(r). This leads to an inhomogeneous Larmor 
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frequency, leading to extraneous dephasing (so-called sec- 
ular relaxation) of transverse spins while it diffuses around 
during the time interval te between the 7r/2 and tt pulses 
in a typical NMR echo measurement. Its effect has been 
studied for a constant gradient (Sen, 1999), a parabolic 
field (i.e. linear field gradient) (poussal,1992) and also 



for a periodic case. ( Bergman, 1995 ). A spherical pore 
and the adjacent spherical dipole source were considered 
by Valckenborg et al ( ,Valckenborg,2003j ). Gillis et al 
( |Gillis,2QQ2l ) noted the strong gradient variation in the 
case of a single spherical source, pointing to more faith- 
ful account of the realistic field profile. Recently, Anand 
et al ( jAnand,2QQ7l ) considered the field produced by an 
arrangement of dipole sources placed on a spherical grain 
and considered qualitatively different regimes. 



gradient derived 




Correlated Noise Texture 





5p/Po 



Figure 14: Two 5p textures imposed on the interface of Figure 



13 From numerically calculated internal field, Bzijo), its lo- 
cal gradient strength \ VBz{y) \ is obtained using the finite dif- 
ference scheme, and SpB is derived from it (shown in the top 
panel along with its probability distribution on the interface). 
The bottom panel shows the texture generated using the corre- 
lated random noise sequence for comparison. In the latter, the 
pattern has a shorter correlation length compared to the former, 
although there are occasional large patch areas of enhanced p 
(e.g. red zones in the bottom corner). The rms-deviations for 
the distributions are 0.65 and 0.35 respectively. 

The inhomogeneous B(r) in real rocks displays some 
aspects not explicitly captured in these studies. For a sim- 
ple geometry, it is obvious that the internal field gradient 
strength depends strongly on the orientation of the inter- 
face. It is further noted that the field variation is strongly 
enhanced closely along the interface with the right ori- 
entation and rapid pore shape variation. (R yu,2Q0 1 ) Start 
with a field profile from an isolated spherical dipole of 
radius a, then extend it for a polydisperse glass beads 
pack as an example of complex pore matrix. From the 
exercise, we hypothesize that the interface is lined with a 
layer in which the innternal field gradient strength is pro- 
nounced with its order of magnitude given by \\/B\ ~ 
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Figure 15: M vs. t for the carbonate pore of Figure 13 for 
different textures of p: uniform po (black), -gradient de- 
rived 6pB (red), and correlated random noise generated texture 
(blue). Short broken lines indicate the predicted initial slope 
with r~'^ = poS/V. The insets show the T2 -distributions with 
the same color scheme. Top panel is for po = 40/im/sec, bot- 
tom panel is for po = 170/im/sec. In both cases, V/S = 



8.88/im, D 
0.23. 



2500/im /sec. Porosity of the sample was 



AxBq I e where e is the typical length scale for the local 
curvature of corrugated pore interface, e may be a small 
fraction of the gain size and depends on the roughness of 
the interface. On a coarser length scale, this rapid varia- 
tion of local field would be smeared out, as the local in- 
terface orientation with respect to the applied field would 
change rapidly. This coarser picture interpolates to what 
one would observe with a dense aggregation of dipoles 
on a smooth spheres. A molecule diffusing in the zone 
over the echo spacing te < ^^/D, however, would feel 
the field variation in a finer scale, and the effect on the 
accumulated phase of the transverse spin, will be equiv- 
alent to that of layer of thickness e lining the part of the 
wall with an surface relaxivity 



PB 



(62) 
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where we assume that jAxBqte < 1 to allow Carr- 
Purcell's classic argument to hold, (otherwise, it crosses 
over to the diffusion controlled regime with weakened 
dependence on te but at values of ps higher than esti- 
mated below.) For = 500G, D = 2500/im/sec^ 
we get 4.46 x iqH ^^^^^ (te in msec, Ax in emu, e 
in /im). For echo spacing of 1 msec, a quartz sand- 
stone with Ax ~ 4 X 10~^ and e = 50/im would give 
Pb ~ 0.1 (/im/s). For matrix with finer length scales 
e < 10/im, and moderate susceptibility contrast Ax ~ 
10~^ (such as realized experimentally via coating of sand 
grains (Keating,2Q07 ) ), pb ~ 100— 1000/im/s is reached. 
Whether this secular mechanism of relaxation overwhelms 
that due to embedded paramagnetic impurities is a ques- 
tion that should be asked for each sample (Ax, e, D) and 
the experimental setup (r^;, Bq), looking for its signature 
in its dependence on varying te- Having shown that it is 
a viable source of an apparent pb within an experimen- 
tally plausible range, we may view the spatial variation 
inherent in AxBo/e as Spb{t^), and ask the same ques- 
tion as to its impact on Aq. 

Detailed analysis of B(r) in a large 3D porous rock is 
beyond the scope of this paper and will be reported else- 
where. Here we briefly summarize how we derive the 
psQudo-Sp texture using the technique we recently devel- 
oped for the internal field and its gradients. The method 
has been tested for a variety of porous media including 
a packed cylinders for which direct experimental imag- 



ing was also done for comparision.(Cho,2009) Figure 12 



shows examples of such calculation applied to a large 3D 
tomogram of carbonate and sandstones. The top panel 
shows the full volumetric rendition of the internal field. 
The bottom panels show cross-section of full 3D calcu- 
lations for the carbonate and the Berea sandstone. To 
minimize the boundary effect, the volume (1.5 x 1.5 x 
1.3cm^) including the surrounding fluid zone was also 
included in the calculation. It was further enlarged by 
mirror reflection images in all three directions. Then the 
mirror-imaged volume was periodically repeated via the 
FFT algorithm used in the calculation of the internal field 
components Ba/{AxBo). The field gradient strength 
was calculated from the finite difference evaluation of 



\2\0.b 



Figure p3] shows the carbonate pore structure used for 
simulating the NMR response using the internal field- 
derived Pb{^) texture as the sole relaxation mechanism. 
To avoid artificial blockage on the bounding surface, a 
sub- volume of 256^ voxels was taken deep inside the 



tomogram( [Sheppard,2004p , then mirrored in all directions. 
Then NMR simulation was performed (Ryu,2008 ) with 
up to 0.8 X 10^ random walkers subject to the stochastic 
dephasing at the wall under the periodic boundary con- 
dition. From the internal field Bz and its gradient, we 



generate the following pseudo- texture: 

PB(r) = Po(l + < |vg^|2 > ) (63) 

where po now corresponds to the effective relaxivity (Eq. 
62 ) arising from a uniform gradient strength with the av- 



erage < \^Bz\ >. The probability distribution is shown 



in the upper panel of Figure 14 For comparison, we also 
applied the correlated-random noise sequence( Ryu,2008| ) 
to generate a patchy Sp texture. Its distribution (lower 
panel of Figure [14]) is symmetric and has about the half 
the rms-deviation as that oiSpB- 

The right panels show how the distribution looks like 
on the pore-grain interface. Note that the | P -derived 
5pB shows a pronounced planar correlation, reflecting 
its sensitivity to the local orientation of the surface with 
respect to the field direction. It also has a skewed distri- 
bution. The carbonate used in this calculation has a pore 
structure somewhat more complex than a sandstone, but 
the grains have relatively smooth surface, therefore we do 
not seem to have strong variation of \VBz\^ on sub-grain 
scales. The correlated random noise texture has varia- 
tions on relatively finer length scales. Given that there is 
a visibly pronounced correlation 'm 5pB on the granular 
length scale, our previous findings would indicate that 
it might have a better chance of having an impact than 
the correlated random-case. They also have contrasting 
skewness in their distributions. In the latter case, unless 
a subtle effect such as the case with quadrature patterns 
considered above, one may expect a gross cancellation 
effect. 

The results of simulation for two different values of 

Note 



15 



po^/D = 0.14,0.60 are shown in Figure 
that in both cases, the curves all have much pronounced 
multi-exponential characteristics compared to the curves 
we have for a simple cubic pore (Figure [TO]). This is un- 
derstandable given the much more complex pore goeme- 
try, yet a naive application of the A>:-criterion would have 
one to expect a more single-exponential like relaxation. 
In all cases, the initial decay faithfully follow the expo- 
nential ex.-p{—tpoS/V) behavior (shown as broken lines) 
as predicted, with and without 5p(r). Note that the range 
of agreement, however, becomes impractically narrow as 
increases. Finally, we observe that the spatial varia- 
tion of p, when it is prescribed to mimic the internal field 
gradient brings about a visible shift in the slowest rate, 
while p with the correlated-random noise pattern hardly 
impacted the relaxation behavior with respect to the uni- 
form case. The T2 distributions, shown in the inset, dis- 
play a hint of minor changes for larger k,, yet it is neg- 
ligible for both cases given the inherent uncertainty in 
the inversion process. For the log interpretation, this im- 
plies that the explanation for the spread and shape of the 
T2 distribution as shown in the figure or as one might 
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have encountered in an NMR logging operation should 
be sought in the geometrical features of the pore matrix 
rather than in Sp{t), as long as there is no compelling rea- 
son to assume an unusually strong magnetic perturbation. 
On the other hand, high quality laboratory data taken with 
controlled variations in the echo spacing, when analyzed 
in the time domain, may display the trend we observe 
numerically. Our findings suggest an interesting possi- 
bility to probe the internal field induced relaxation us- 
ing artificial bead/grain packs in which one can vary the 
smoothness and anisotropy of individual grains and their 
arrangement. 

CONCLUSIONS 

Using the rectangular pore geometry, we developed an 
analytical method to systematically probe the intertwined 
effect of its geometrical and p textural variations. The 
result yields the fractional changes in the slowest relax- 
ation rate with respect to that of the uniform p case for 
various values of and a which control the geometry 
and p— texture respectively. This framework provides a 
bound for the uncertainty in the NMR log interpretation 
for complex formations. We identified the relative sym- 
metry of the slowest eigenmode and the interfacial pat- 
tern of as factors that control the impact of varia- 
tion, and demonstrated the mechanism using a face cen- 
tered and corner centered Spin a. cubic pore. While it is 
plausible that the deposit pattern of paramagnetic ions 
may be incommensurate with the underlying pore ge- 
ometry, it is not unthinkable that some correlation exist 
due to the diagenesis process constrained by the complex 
pore geometry. Currently, there is not empirical data with 
enough details to further narrow down on the p texture. 
Experimental systems using coated glass bead packs may 
be designed to shed further light on the issue. Colloidal 
suspension with controllable interfacial properties is an- 
other system where our finding may prove useful. In real 
rocks, detailed mapping of magnetic profile of the pore- 
grain interface on small enough length scales would be 
invaluable for further progress. 

In the absence of such detailed empirical constraint in 
a typical NMR logging operation, we use as guidance 
simulation results based on artificially made 6p textures 
imposed on a tomogram-derived pore structure. Prevously, 
we demonstrated that the salient features observed on 
the simple closed-pore or bead packs remain in tact for 
a sandstone. There, we made comparisons among tex- 
tures of Sp, quasi randomly generated, with different de- 
grees of spatial correlation. In the current work, we fur- 
ther made comparisons of a different nature, between one 
generated from the internal field gradient, and the other 
by the random noise sequence. The result suggests that 
gross cancelation due to symmetric grounds or short cor- 



relation lengths of the texture seem to mitigate the impact 
of (^p in a typical natural media, while stronger spatial 
correlation as in the internal field gradient profile, tends 
to act against such a trend leading to a weak, but observ- 
able shift. 

The phase- space is too vast for drawing any conclusive 
generalization out of these numerical results alone. Still, 
in view of analytical results we elaborated in this work, 
it seems reasonable to conclude that the effect of Sp on 
logging application is secondary only to the more signif- 
icant effect of pore geometry variation as far as the dom- 
inant features of the T2 distribution are concerned. The 
fractional shift and its inferred bounds (Figs|4]|6]), when 
translated into uncertainties in the pore length scale, pro- 
vide useful guidance when there is not enough informa- 
tion to determine the relative significance of the geomet- 
rical and the lithological variation in a given formation. 
While its impact seems minimal in the overall T2 distri- 
bution of the NMR logging data, it is worth pointing out 
that the transport properties tend to be controlled by fea- 
tures of the matrix which are manifest only in a narrow 
range of scales. The effect of internal field gradient, for 
having close relationship with pore geometry variations 
and also for inducing visible changes in NMR response, 
may be utilized toward such purpose once we gain better 
understanding and control over its impact. Therefore, in 
our effort to refine the NMR logging as a more effective 
permeability probe, detailed aspects of Sp and its impact 
on NMR response should merit further investigations. 
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